Preliminary Bayesian Parameter Estimation Model

TODO: write an introductory section TODO: add commentary throughout to guide reader through the steps of the tutorial

# set filepath for data file
filepath <- "https://raw.githubusercontent.com/LRI-2/Data/main/ILD/AMIBshare_persons_2019_0501.csv"
# read in the .csv file using the url() function
AMIB_persons <- read.csv(file=url(filepath),header=TRUE)
# set filepath for data file
filepath <- "https://raw.githubusercontent.com/LRI-2/Data/main/ILD/AMIBshare_daily_2019_0501.csv"
# read in the .csv file using the url() function
AMIB_daily <- read.csv(file=url(filepath),header=TRUE)
# set filepath for data file
filepath <- "https://raw.githubusercontent.com/LRI-2/Data/main/ILD/AMIBshare_interaction_2019_0501.csv"
# read in the .csv file using the url() function
AMIB_interaction <- read.csv(file=url(filepath),header=TRUE)
# subset to day-level variables of interest
bpe_daily <- AMIB_daily[,c("id", "day", "slphrs", "negaff")]
# subset to person-level variables of interest
bpe_persons <- AMIB_persons[,c("id", "erq_reap", "erq_supp")]
# merge day- and person-level data
bpe_data <- merge(bpe_daily, bpe_persons, by = "id")
# center the erq subscale score variables for interpretability
bpe_data$erq_reap_c <- scale(bpe_data$erq_reap, center=TRUE,scale=FALSE)[,1]
bpe_data$erq_supp_c <- scale(bpe_data$erq_supp, center=TRUE,scale=FALSE)[,1]
#TODO: plot the raw data to get a sense for any trends in the data
      # correlation matrix
      # person-by-person plot
bpe.reap <- brm(bf(negaff ~ b0 + b1 * day + b2 * slphrs + b3 * erq_reap_c + b4 * slphrs * erq_reap_c, b0 ~ 1 + (1|id), b1 ~ 1 + (1|id), b2 ~ 1 + (1|id), b3 ~ 1 + (1|id), b4 ~ 1 + (1|id), nl = TRUE), data = bpe_data, family = gaussian(), prior = c(
  prior(normal(0, 1), nlpar = "b0"),
  prior(normal(0, 1), nlpar = "b1"),
  prior(normal(0, 1), nlpar = "b2"),
  prior(normal(0, 1), nlpar = "b3"),
  prior(normal(0, 1), nlpar = "b4")),
  iter = 2000, chains = 4, cores = 4)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## Warning: There were 2 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
summary(bpe.reap)
## Warning: There were 2 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: negaff ~ b0 + b1 * day + b2 * slphrs + b3 * erq_reap_c + b4 * slphrs * erq_reap_c 
##          b0 ~ 1 + (1 | id)
##          b1 ~ 1 + (1 | id)
##          b2 ~ 1 + (1 | id)
##          b3 ~ 1 + (1 | id)
##          b4 ~ 1 + (1 | id)
##    Data: bpe_data (Number of observations: 1421) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~id (Number of levels: 190) 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(b0_Intercept)     0.58      0.05     0.47     0.69 1.00      882     1786
## sd(b1_Intercept)     0.07      0.02     0.03     0.09 1.00      597      398
## sd(b2_Intercept)     0.03      0.01     0.00     0.05 1.02      298      750
## sd(b3_Intercept)     0.10      0.07     0.00     0.27 1.00      822     1330
## sd(b4_Intercept)     0.01      0.01     0.00     0.04 1.00      931     1344
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## b0_Intercept     3.21      0.11     2.99     3.42 1.00     4089     2986
## b1_Intercept    -0.07      0.01    -0.09    -0.05 1.00     5836     3208
## b2_Intercept    -0.07      0.01    -0.10    -0.05 1.00     4020     2736
## b3_Intercept    -0.02      0.10    -0.22     0.18 1.00     4128     3054
## b4_Intercept     0.00      0.01    -0.03     0.03 1.00     4817     3064
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.77      0.02     0.74     0.81 1.00     1963     2716
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayes_R2(bpe.reap)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.4535433 0.01942911 0.4149674 0.4898165
plot(bpe.reap)

TODO: interpret model results (including Bayesian significance testing) TODO: make a visualization of this model that shows the prototypical individual and all individuals

bpe.supp <- brm(bf(negaff ~ b0 + b1 * day + b2 * slphrs + b3 * erq_supp_c + b4 * slphrs * erq_supp_c, b0 ~ 1 + (1|id), b1 ~ 1 + (1|id), b2 ~ 1 + (1|id), b3 ~ 1 + (1|id), b4 ~ 1 + (1|id), nl = TRUE), data = bpe_data, family = gaussian(), prior = c(
  prior(normal(0, 1), nlpar = "b0"),
  prior(normal(0, 1), nlpar = "b1"),
  prior(normal(0, 1), nlpar = "b2"),
  prior(normal(0, 1), nlpar = "b3"),
  prior(normal(0, 1), nlpar = "b4")),
  iter = 2000, chains = 4, cores = 4)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.0.40.1)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
## Warning: There were 6 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
summary(bpe.supp)
## Warning: There were 6 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: negaff ~ b0 + b1 * day + b2 * slphrs + b3 * erq_supp_c + b4 * slphrs * erq_supp_c 
##          b0 ~ 1 + (1 | id)
##          b1 ~ 1 + (1 | id)
##          b2 ~ 1 + (1 | id)
##          b3 ~ 1 + (1 | id)
##          b4 ~ 1 + (1 | id)
##    Data: bpe_data (Number of observations: 1421) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~id (Number of levels: 190) 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(b0_Intercept)     0.57      0.06     0.45     0.68 1.00      653     1367
## sd(b1_Intercept)     0.06      0.01     0.04     0.09 1.00      896      994
## sd(b2_Intercept)     0.03      0.01     0.00     0.06 1.01      266      726
## sd(b3_Intercept)     0.11      0.07     0.00     0.26 1.01      522     1167
## sd(b4_Intercept)     0.01      0.01     0.00     0.03 1.00      702      945
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## b0_Intercept     3.19      0.12     2.96     3.42 1.00     2430     2249
## b1_Intercept    -0.07      0.01    -0.09    -0.05 1.00     6258     3312
## b2_Intercept    -0.07      0.01    -0.10    -0.04 1.00     2569     2650
## b3_Intercept     0.12      0.08    -0.04     0.29 1.00     3650     2984
## b4_Intercept    -0.01      0.01    -0.03     0.01 1.00     4271     3114
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.77      0.02     0.74     0.81 1.00     2308     2708
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayes_R2(bpe.supp)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.4536001 0.01878353 0.4156904 0.4887272
plot(bpe.supp)

TODO: interpret model results (including Bayesian significance testing) TODO: make a visualization of this model that shows the prototypical individual and all individuals

TODO: compare the two models (reap vs. supp - which one seems more predictive of better control of negaff?)

TODO: write up a final analysis/conclusion for this tutorial